Put all scenarios in a loop

This notebook is used to visualize the congested lines in WECC. It is adapted from L:\Renewable Energy\EnergyGridModeling\InteractiveInterface\InteractiveMapVersion01.ipynb

TO DO:

  • use the zscore that was calculated in the congestion_pvalues notebook: input congestion file should include z_score and p_value columns (need both: zscore is just for normal approx, pvalues work for any dist)
  • zip colors and linewidths together
  • include the full WECC list of regions
  • map other details (renewable energy sources, etc.)
  • once WECC data are available, include ability to highlight differences
In [1]:
import numpy as np
import pandas as pd

import matplotlib.cm as cm
import matplotlib

from bokeh.io import show,output_notebook,output_file
from bokeh.models import (
    ColumnDataSource,
    HoverTool,
    LogColorMapper,
    CustomJS,
    Slider,
    CheckboxGroup,
    Select,
    Div,
    TapTool
)
from bokeh.palettes import Viridis6 as palette
from bokeh.plotting import figure
from bokeh.layouts import (widgetbox, row)
from bokeh.tile_providers import CARTODBPOSITRON,STAMEN_TONER_BACKGROUND,STAMEN_TONER
from bokeh.sampledata.us_counties import data as counties
from bokeh import events
from bokeh.colors.groups import yellow

import westernintnet
In [2]:
wecc = westernintnet.WesternIntNet()
Loading sub
Loading bus2sub
Loading bus
Loading genbus
Loading branches
Loading resources
Loading net_generation
Load solar data
Load wind data
Load hydro data
Load demand data
Done loading
In [3]:
output_notebook()
Loading BokehJS ...

Projection from lon/lat to marcator for openstreetmap

In [4]:
from pyproj import Proj, transform

def reproject_wgs_to_itm(x_lon_lat):
    prj_wgs = Proj(init='epsg:4326')
    prj_itm = Proj(init='EPSG:3857')
    x, y = transform(prj_wgs, prj_itm, x_lon_lat[0], x_lon_lat[1])
    r = [x,y]
    return r

County names and shapes

In [5]:
state_list = ['wa','or','ca','nv','az','ut','nm','co','wy','id','mt','tx']

#This just shades CA for reference; 
counties = {
    code: county for code, county in counties.items() if county["state"] == "ca"
}

#This ends up giving the whole world
#counties = {
#    code: county for code, county in counties.items() for county["state"] in state_list
#}

# Convert lon/lat to x/y
county_xys = [reproject_wgs_to_itm([county["lons"],county["lats"]]) for county in counties.values()]

county_xs = [county_xys[i][0] for i,county in enumerate(counties.values())]
county_ys = [county_xys[i][1] for i,county in enumerate(counties.values())]

county_names = [county['name'] for county in counties.values()]
county_color = []
for c_name in county_names:
    color = 'black'
    county_color.append(color) 

source = ColumnDataSource(data=dict(
    x=county_xs,
    y=county_ys,
    name=county_names,
    color=county_color,
))
In [6]:
r1 = wecc.branches[['from_lon','from_lat']].apply(reproject_wgs_to_itm,axis=1)
wecc.branches['from_x'] = r1.apply(lambda x:x[0])
wecc.branches['from_y'] = r1.apply(lambda x:x[1])

r2 = wecc.branches[['to_lon','to_lat']].apply(reproject_wgs_to_itm,axis=1)
wecc.branches['to_x'] = r2.apply(lambda x:x[0])
wecc.branches['to_y'] = r2.apply(lambda x:x[1])

wecc_base_branches = wecc.branches.copy()
In [7]:
#Add steps to fix branch capacity for the congested branches 2016,2020,2030
In [8]:
def makeMap(congestion_df):
    util_rate = congestion_df['hutil>=0p75']
    wecc_branches  = pd.concat([wecc_base_branches, congestion_df], axis=1)

    #Replace this population-derived values from congestion_p-values notebook, or use the zsore from the notebook remove this
    mean = wecc_branches.loc[wecc_branches['Capacity']!=99999]['hutil>=0p75'].mean()
    sdev = wecc_branches.loc[wecc_branches['Capacity']!=99999]['hutil>=0p75'].std()

    minima = util_rate.min()
    maxima = util_rate.max()

    norm =  matplotlib.colors.Normalize(vmin=minima, vmax=maxima, clip=True)
    mapper = cm.ScalarMappable(norm=norm, cmap=cm.jet)
    mapper.set_array([])
    colors = mapper.to_rgba(util_rate.values)

    # Pick colors and linewidths for 3 types of lines depending on zscore = (x-mu)/sigma
    # gray, light red, dark red for increasing congestion
    colors = ['#AAB7BB','#FF5733','#C70039']  
    # thin, thick, thick for increasing congestion
    linewidths = [0.5,5.0,5.0]  

    #replacing with pre-calculated z-score
    zscore_t1 = 2.
    zscore_t2 = 1.
    wecc_branches.loc[(wecc_branches['zscore'] >= zscore_t1),'color'] = colors[2]
    wecc_branches.loc[(wecc_branches['zscore']  < zscore_t1) 
                      & (wecc_branches['zscore']  >= zscore_t2),'color'] = colors[1]
    wecc_branches.loc[(wecc_branches['zscore'] < zscore_t2),'color'] = colors[0]

    wecc_branches.loc[(wecc_branches['zscore']  >= zscore_t1),'linewidth'] = linewidths[2]
    wecc_branches.loc[(wecc_branches['zscore']  < zscore_t1) 
                      & (wecc_branches['zscore']  >= zscore_t2),'linewidth'] = linewidths[1]
    wecc_branches.loc[(wecc_branches['zscore']  < zscore_t2),'linewidth'] = linewidths[0]

    #set up figure
    TOOLS = "pan,wheel_zoom,reset,hover,save"

    p = figure(
        title="WECC Energy Grid", tools=TOOLS,
        x_axis_location=None, y_axis_location=None,
        plot_width=800,plot_height=800)

    p.add_tile(CARTODBPOSITRON)
    p.patches('x', 'y', source=source, fill_color={'field':'color'}, fill_alpha=0.1, line_color="white",  line_width=0.5)
    p.multi_line(wecc_branches[['from_x','to_x']].values.tolist(), 
                 wecc_branches[['from_y','to_y']].values.tolist(), 
                 line_color=wecc_branches['color'],line_width=wecc_branches['linewidth'])
    #show(p)
    return p
In [9]:
datadir = 'U:\\src\\PostREISE\\postreise\\analyze\\data\\'
scenarios = ['western_scenario_Update01','california2020Test01','california2030Test01',
          'california2020_fixCalCong','california2030_fixCalCong','california2020_westTarget','california2030_westTarget',
           'california2020_fixCalCongX','california2030_fixCalCongX' ]

#scenarios = ['western_scenario_Update01','california2020Test01','california2030Test01',
#          'california2020_westTarget','california2030_westTarget']

congestion_map = {}
for s in scenarios:
    congestion_source = datadir + 'congestion_' + s +'.csv'
    congestion_df = pd.read_csv(congestion_source)
    congestion_map[s] = makeMap(congestion_df)

for s in scenarios:
    print(s)
    show(congestion_map[s])
western_scenario_Update01
california2020Test01
california2030Test01
california2020_fixCalCong
california2030_fixCalCong
california2020_westTarget
california2030_westTarget
california2020_fixCalCongX
california2030_fixCalCongX
In [10]:
wecc.branches.head()
Out[10]:
fbus tbus r x b rateA rateB rateC ratio angle ... mu_angmin mu_angmax from_lat from_lon to_lat to_lon from_x from_y to_x to_y
0 10002 10001 0.057474 0.297874 0.04557 185.33 0 0 0.0 0.0 ... 0.0 0.0000; 47.695555 -124.183645 48.241400 -124.577777 -1.382406e+07 6.056355e+06 -1.386793e+07 6.147110e+06
1 10011 10001 0.027895 0.200986 0.06271 167.88 0 0 0.0 0.0 ... 0.0 0.0000; 48.002536 -123.762001 48.241400 -124.577777 -1.377712e+07 6.107277e+06 -1.386793e+07 6.147110e+06
2 10014 10002 0.028286 0.132954 0.07205 207.43 0 0 0.0 0.0 ... 0.0 0.0000; 47.188929 -123.686007 47.695555 -124.183645 -1.376866e+07 5.972967e+06 -1.382406e+07 6.056355e+06
3 10004 10003 0.012099 0.073288 0.01171 170.67 0 0 0.0 0.0 ... 0.0 0.0000; 46.927528 -124.171950 47.040041 -124.056969 -1.382276e+07 5.930253e+06 -1.380996e+07 5.948612e+06
4 10003 10010 0.012423 0.079127 0.01766 201.54 0 0 0.0 0.0 ... 0.0 0.0000; 47.040041 -124.056969 47.174349 -123.847397 -1.380996e+07 5.948612e+06 -1.378663e+07 5.970579e+06

5 rows × 29 columns

In [11]:
wecc.branches.columns
Out[11]:
Index(['fbus', 'tbus', 'r', 'x', 'b', 'rateA', 'rateB', 'rateC', 'ratio',
       'angle', 'status', 'angmin', 'angmax', 'Pf', 'Qf', 'Pt', 'Qt', 'mu_Sf',
       'mu_St', 'mu_angmin', 'mu_angmax', 'from_lat', 'from_lon', 'to_lat',
       'to_lon', 'from_x', 'from_y', 'to_x', 'to_y'],
      dtype='object')